Mixed Effects Models: Effect of Soil Moisture Variables
Soil moisture is a covariate that may impact seed production. In the field, we measured soil moisture in designated areas of the plot once per week, for 9 weeks. To understand which variables representing change in soil moisture are most influential for seed production, we ran a generalized linear mixed model with all of the soil moisture variables. Site, plot treatment, and species are random effects, while the soil moisture variables (mean moisture, rate of change, effective minimum or value at week 7, final moisture reading, and range) are fixed effects. Once we had this model, we ran the dredge function to find the most appropriate combination of variables for our larger model.
We are including all of the individuals of the three species and using species as a random effect.
First, we combined the variables for the three species in order to run the model selection with all of the collected seeds. We checked the distribution of the total number of developed seeds and found that they fit a negative binomial distribution. Next we checked the distribution for seeds per flower and found that they fit an exponential distribution.
condensed.mert<-mert[,c(-4,-20,-21,-23,-24,-30)] #removing unnecessary columns
condensed.delph<-delph[,c(-4,-20:-28,-30,-36)]
condensed.pot<-pot[,c(-4,-20:-28,-30,-33:-37)]
names(condensed.mert)[names(condensed.mert)=="dev.seed"]<-"total.seeds"
condensed.mert<-condensed.mert[(condensed.mert$treat=="open"),] #limiting to open treatment
condensed.delph<-condensed.delph[(condensed.delph$treat=="open"),]
condensed.pot<-condensed.pot[(condensed.pot$treat=="open"),]
all.seeds<-rbind(condensed.mert,condensed.delph,condensed.pot) #combining all of the data sets
all.seed.nbinom<-fitdist(all.seeds$total.seeds, "nbinom") #fitting Mertensia data to negative binomial distribution using fitdistplus
plot(all.seed.nbinom) #plotting to see the fit
#all.seeds.exp<-fitdist(all.seeds$seeds.per.flower,"exp") #checking exponential distribution
#plot(all.seeds.exp) #plotting exponential distribution
Because the total seeds fit a negative binomial distribution, we used the glmmTMB package. Because the seeds per flower variable fit a exponential distribution, we used the lme4 package. The seeds per flower variable requires the Gamma family, which means we had to create a hurdle model to avoid the zeros in our data. Hurdle models involve a first binomial model with zeros and non-zeros as the response variable, and a second Gamma model with seeds per flower as the response variable. The first model shows us if the soil moisture impacts the ability for species to produce seeds (success/failure).
soil.mod<-glmmTMB(total.seeds ~ mean + rate + effect.min + end + range + (1|site/plot.treat) + (1|species), family = "nbinom2", data = all.seeds) #model for soil moisture variables with all seeds
summary(soil.mod)
## Family: nbinom2 ( log )
## Formula:
## total.seeds ~ mean + rate + effect.min + end + range + (1 | site/plot.treat) +
## (1 | species)
## Data: all.seeds
##
## AIC BIC logLik deviance df.resid
## 1331.8 1360.5 -655.9 1311.8 120
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plot.treat:site (Intercept) 1.387e-08 1.178e-04
## site (Intercept) 2.567e-11 5.066e-06
## species (Intercept) 5.321e-01 7.295e-01
## Number of obs: 130, groups: plot.treat:site, 14; site, 7; species, 3
##
## Overdispersion parameter for nbinom2 family (): 1.29
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.97353 0.59194 5.023 5.08e-07 ***
## mean -0.05567 0.06613 -0.842 0.3999
## rate 0.03893 0.13544 0.287 0.7738
## effect.min 0.13025 0.07335 1.776 0.0758 .
## end 0.02627 0.03829 0.686 0.4927
## range 0.07068 0.07024 1.006 0.3143
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
all.seeds$non_zero <- ifelse(all.seeds$seeds.per.flower > 0, 1, 0) # creating zero/non-zero column
soil.mod2.hurdle<-glmer(non_zero ~ mean + rate + effect.min + end + range + (1|site/plot.treat) + (1|species), data=all.seeds,family = binomial) # model with new column
all.seeds.nonzero = subset(all.seeds,non_zero==1)
soil.mod2<-glmer(seeds.per.flower ~ mean + rate + effect.min + end + range + (1|site/plot.treat) + (1|species), family=Gamma,data = all.seeds.nonzero,na.action = "na.fail") # model without zeros
summary(soil.mod2.hurdle) # viewing model results
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula:
## non_zero ~ mean + rate + effect.min + end + range + (1 | site/plot.treat) +
## (1 | species)
## Data: all.seeds
##
## AIC BIC logLik deviance df.resid
## 69.0 94.0 -25.5 51.0 109
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -5.7889 0.1951 0.2221 0.2667 0.4469
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot.treat:site (Intercept) 0 0
## site (Intercept) 0 0
## species (Intercept) 0 0
## Number of obs: 118, groups: plot.treat:site, 14; site, 7; species, 3
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 1.9004758 2.0806674 0.913 0.361
## mean 0.0345014 0.3093144 0.112 0.911
## rate -0.0006426 0.6011741 -0.001 0.999
## effect.min 0.0091411 0.3367457 0.027 0.978
## end 0.1447328 0.1770694 0.817 0.414
## range -0.1040140 0.3284441 -0.317 0.751
##
## Correlation of Fixed Effects:
## (Intr) mean rate effct. end
## mean -0.196
## rate -0.038 -0.859
## effect.min -0.161 -0.798 0.725
## end -0.400 -0.074 0.122 -0.166
## range -0.002 -0.928 0.950 0.795 0.016
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
summary(soil.mod2)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( inverse )
## Formula: seeds.per.flower ~ mean + rate + effect.min + end + range + (1 |
## site/plot.treat) + (1 | species)
## Data: all.seeds.nonzero
##
## AIC BIC logLik deviance df.resid
## 525.0 552.1 -252.5 505.0 101
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.6365 -0.6956 -0.1135 0.5339 2.7921
##
## Random effects:
## Groups Name Variance Std.Dev.
## plot.treat:site (Intercept) 2.224e-09 4.716e-05
## site (Intercept) 1.211e-10 1.101e-05
## species (Intercept) 1.028e-01 3.206e-01
## Residual 3.464e-01 5.885e-01
## Number of obs: 111, groups: plot.treat:site, 14; site, 7; species, 3
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 0.5528490 0.3166158 1.746 0.0808 .
## mean 0.0036499 0.0039235 0.930 0.3522
## rate 0.0100715 0.0075076 1.341 0.1798
## effect.min -0.0086623 0.0056242 -1.540 0.1235
## end -0.0006279 0.0034459 -0.182 0.8554
## range -0.0016415 0.0035777 -0.459 0.6464
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) mean rate effct. end
## mean -0.044
## rate -0.022 -0.450
## effect.min -0.018 -0.628 0.478
## end -0.035 -0.019 -0.077 -0.416
## range -0.002 -0.685 0.874 0.570 -0.169
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
In our soil moisture hurdle model, the soil moisture variables did not affect the success or failure of the plants to produce seeds (zero-nonzero), and a few soil moisture variables did affect the number of seeds/flower. Therefore, we can run a model selection on the second model that excludes the zeros.
We ran two models selections: one for the glmmTMB model with the count data and one for the glmer model with the seeds/flower data.
dredge_result<-dredge(soil.mod) #using dredge function for model selection
subset(dredge_result, delta<4)
## Global model call: glmmTMB(formula = total.seeds ~ mean + rate + effect.min + end +
## range + (1 | site/plot.treat) + (1 | species), data = all.seeds,
## family = "nbinom2", ziformula = ~0, dispformula = ~1)
## ---
## Model selection table
## cnd((Int)) dsp((Int)) cnd(eff.min) cnd(end) cnd(men) cnd(rng) cnd(rat) df
## 18 3.248 + 0.07976 -0.09793 7
## 10 3.020 + 0.08948 0.033050 7
## 17 3.892 + -0.10040 6
## 14 3.104 + 0.12880 -0.040740 0.055870 8
## 2 3.558 + 0.06578 6
## 21 3.447 + 0.028800 -0.09829 7
## 3 3.507 + 0.05534 6
## 1 4.084 + 5
## 9 3.757 + 0.032290 6
## 26 3.064 + 0.08903 0.015690 -0.06937 8
## 19 3.530 + 0.03764 -0.08212 7
## 20 3.126 + 0.07006 0.02068 -0.08921 8
## 22 3.149 + 0.07663 0.007909 -0.09739 8
## 12 2.867 + 0.07883 0.02558 0.030120 8
## 5 3.589 + 0.031530 6
## 4 3.256 + 0.04996 0.04104 7
## 25 3.833 + 0.008667 -0.08560 7
## 11 3.383 + 0.04393 0.023500 7
## 16 2.996 + 0.11410 0.02314 -0.038910 0.051420 9
## 6 3.368 + 0.05556 0.017280 7
## 7 3.319 + 0.04660 0.017700 7
## 29 3.365 + 0.045060 -0.026080 -0.14530 8
## 23 3.302 + 0.02452 0.022490 -0.09045 8
## 30 3.102 + 0.13450 -0.046090 0.062150 0.01232 9
## 13 3.598 + 0.014170 0.026020 7
## 28 2.960 + 0.07902 0.01960 0.015020 -0.06181 9
## logLik AICc delta weight
## 18 -656.608 1328.1 0.00 0.107
## 10 -656.910 1328.7 0.60 0.079
## 17 -658.035 1328.8 0.62 0.079
## 14 -656.162 1329.5 1.38 0.054
## 2 -658.492 1329.7 1.53 0.050
## 21 -657.430 1329.8 1.65 0.047
## 3 -658.556 1329.8 1.66 0.047
## 1 -659.679 1329.8 1.71 0.046
## 9 -658.639 1330.0 1.83 0.043
## 26 -656.409 1330.0 1.88 0.042
## 19 -657.550 1330.0 1.89 0.042
## 20 -656.461 1330.1 1.98 0.040
## 22 -656.551 1330.3 2.16 0.036
## 12 -656.690 1330.6 2.44 0.032
## 5 -658.968 1330.6 2.49 0.031
## 4 -657.923 1330.8 2.63 0.029
## 25 -657.994 1330.9 2.77 0.027
## 11 -658.004 1330.9 2.79 0.027
## 16 -655.965 1331.4 3.30 0.021
## 6 -658.277 1331.5 3.34 0.020
## 7 -658.330 1331.6 3.44 0.019
## 29 -657.233 1331.7 3.52 0.018
## 23 -657.249 1331.7 3.55 0.018
## 30 -656.158 1331.8 3.68 0.017
## 13 -658.540 1332.0 3.87 0.016
## 28 -656.275 1332.0 3.92 0.015
## Models ranked by AICc(x)
## Random terms (all models):
## 'cond(1 | site/plot.treat)', 'cond(1 | species)'
dredge_result2<-dredge(soil.mod2) #using dredge function for model selection
subset(dredge_result2, delta<4)
## Global model call: glmer(formula = seeds.per.flower ~ mean + rate + effect.min +
## end + range + (1 | site/plot.treat) + (1 | species), data = all.seeds.nonzero,
## family = Gamma, na.action = "na.fail")
## ---
## Model selection table
## (Int) eff.min end men rng rat df logLik AICc
## 18 0.5682 -0.005889 0.012260 7 -252.859 520.8
## 17 0.5237 0.013150 6 -254.040 520.9
## 22 0.5501 -0.007555 0.0021860 0.013470 8 -252.562 522.5
## 19 0.5483 -0.0026230 0.012740 7 -253.757 522.6
## 25 0.5110 0.002209 0.018610 7 -253.788 522.7
## 26 0.5571 -0.005505 0.001487 0.015870 8 -252.731 522.9
## 14 0.5541 -0.012990 0.0070630 -0.006219 8 -252.808 523.0
## 20 0.5693 -0.005718 -0.0002574 0.012200 8 -252.856 523.1
## 21 0.5297 -0.0004193 0.012810 7 -254.026 523.1
## 10 0.5914 -0.007620 -0.002983 7 -254.406 523.9
## 27 0.5360 -0.0029340 0.002831 0.019200 8 -253.279 524.0
## 1 0.5028 5 -256.726 524.0
## 2 0.5585 -0.007074 6 -255.637 524.1
## 29 0.5371 -0.0027120 0.004387 0.021820 8 -253.436 524.3
## 9 0.5297 -0.002745 6 -255.780 524.4
## 3 0.5420 -0.0038140 6 -255.984 524.8
## 30 0.5507 -0.009253 0.0037950 -0.001933 0.009577 9 -252.510 524.8
## delta weight
## 18 0.00 0.169
## 17 0.08 0.162
## 22 1.73 0.071
## 19 1.80 0.069
## 25 1.86 0.067
## 26 2.07 0.060
## 14 2.22 0.056
## 20 2.32 0.053
## 21 2.33 0.053
## 10 3.09 0.036
## 27 3.17 0.035
## 1 3.22 0.034
## 2 3.28 0.033
## 29 3.48 0.030
## 9 3.56 0.028
## 3 3.97 0.023
## 30 4.00 0.023
## Models ranked by AICc(x)
## Random terms (all models):
## '1 | site/plot.treat', '1 | species'
After running the model selections, we found that rate of soil moisture with effective minimum moisture was the best model for the seeds/flower response variable and the end soil moisture was the best model for the seed count data.
Mixed Effects Models with Count Data
Generalized linear mixed effects models were appropriate for this analysis because the data contain sources of non-independence. Individual plants were located in the same site as other individuals and were likely genetically similar. Therefore, we considered site a random effect.
The response variable is total seeds, which is discrete and non-zero. The interaction between deviation and the early/late factor was a fixed effect, as well as the interaction between plot treatment and conspecific density. We would have included soil moisture as a covariate and fixed effect in our model, but soil moisture was not important for seed production. Each species is analyzed separately.
Mertensia
We started by limiting the individuals to the open treatment (excluding the hand-pollinated treatment) and by checking the distribution of the total seed counts.
mert.open<-mert[(mert$treat=="open"),] #removed hand pollination treatment
mert.nbinom<-fitdist(mert.open$dev.seed, "nbinom") #fitting Mertensia data to negative binomial distribution using fitdistplus
plot(mert.nbinom) #plotting to see the fit
The Mertensia seed data fits a negative binomial distribution. The glmmTMB package is a good package for modeling data with a negative binomial. The glmmTMB package is also suitable for data that could be overdispersed or zero-inflated.
Main Model
#mert.glmmtmb<-glmmTMB(dev.seed ~ flowers:(deviation*early.late) + flowers:number.conspecifics/plot.treat + (1|site), family = "nbinom2", data = mert.open) #putting count data into glmmTMB for DHARMa package
#summary(mert.glmmtmb) #summary of glmmTMB model
#mert.counts.simulation<-simulateResiduals(fittedModel = mert.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
#plot(mert.counts.simulation) #plotting simulation
#testDispersion(mert.counts.simulation) #checking for overdispersion
#testZeroInflation(mert.counts.simulation) #checking for zero inflation
We detected a significant effect of early vs. late blooming on the total number of seeds in Mertensia individuals, but when we checked for overdispersion and zero-inflation in the data using the DHARMa package, the data were slightly zero-inflated. To address this, we adjusted our glmmtmb model to include zero-inflation.
mert.glmmtmb<-glmmTMB(dev.seed ~ flowers:(deviation*early.late) + flowers:number.conspecifics/plot.treat + (1|site), ziformula=~1, family = "nbinom2", data = mert.open) #putting count data into glmmTMB for DHARMa package
summary(mert.glmmtmb) #summary of glmmTMB model
## Family: nbinom2 ( log )
## Formula:
## dev.seed ~ flowers:(deviation * early.late) + flowers:number.conspecifics/plot.treat +
## (1 | site)
## Zero inflation: ~1
## Data: mert.open
##
## AIC BIC logLik deviance df.resid
## 386.2 402.4 -184.1 368.2 36
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.07728 0.278
## Number of obs: 45, groups: site, 4
##
## Overdispersion parameter for nbinom2 family (): 3.03
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) 2.613e+00 3.761e-01 6.949
## flowers:deviation 2.433e-03 1.540e-03 1.579
## flowers:early.latelate 1.931e-02 1.009e-02 1.914
## flowers:number.conspecifics -2.368e-05 8.772e-05 -0.270
## flowers:deviation:early.latelate -2.693e-03 2.934e-03 -0.918
## flowers:number.conspecifics:plot.treatmanip 2.973e-04 1.937e-04 1.535
## Pr(>|z|)
## (Intercept) 3.68e-12 ***
## flowers:deviation 0.1143
## flowers:early.latelate 0.0557 .
## flowers:number.conspecifics 0.7872
## flowers:deviation:early.latelate 0.3587
## flowers:number.conspecifics:plot.treatmanip 0.1248
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.3800 0.5513 -4.317 1.58e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mert.counts.simulation<-simulateResiduals(fittedModel = mert.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
plot(mert.counts.simulation) #plotting simulation
testDispersion(mert.counts.simulation) #checking for overdispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.90351, p-value = 0.784
## alternative hypothesis: two.sided
testZeroInflation(mert.counts.simulation) #checking for zero inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0695, p-value = 1
## alternative hypothesis: two.sided
When number of flowers is taken into account, seed set for Mertensia individuals may be affected by blooming before or after the population peak. Zero-inflation is accounted for in the new model, and the model for Mertensia count data is no longer zero-inflated.
To figure out the direction of the effect, we looked at the plots (shown below). Blooming after the population peak may increase seed set for Mertensia individuals.
Pollen Limitation Model
To test the pollen limitation individuals, we created a separate fixed effects model. This model has site and plot treatment as random effects and the interaction between pollen limitation treatment and number of flowers as a fixed effect. We included this interaction because the difference in seed set between the open and hand treatments is heavily dependent on the number of flowers. The interaction accounts for the wide variation of number of flowers and ultimately developed seeds. The response variable is still total number of developed seeds in Mertensia individuals.
#mert.treat.glmmtmb<-glmmTMB(dev.seed ~ treat/flowers + (1|site/plot.treat), family = "nbinom2", data = mert) #Mertensia model for negative binomial data with glmmTMB package
#summary(mert.treat.glmmtmb) #summary of model output
#mert.counts.treat.simulation<-simulateResiduals(fittedModel = mert.treat.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
#plot(mert.counts.treat.simulation) #plotting simulation
#testDispersion(mert.counts.treat.simulation) #checking for overdispersion
#testZeroInflation(mert.counts.treat.simulation) #checking for zero inflation
We detected an effect of pollen limitation treatment on Mertensia total developed seeds, but we had to test the model assumptions of our pollen limitation model as well. The data in the pollen limitation model for Mertensia are not overdispersed, but they are zero-inflated. We then adjusted our model to include zero-inflation and ran the overdispersion and zero-inflation tests again.
mert.treat.glmmtmb<-glmmTMB(dev.seed ~ treat/flowers + (1|site/plot.treat), ziformula=~1, family = "nbinom2", data = mert) #Mertensia model for negative binomial data with glmmTMB package
summary(mert.treat.glmmtmb) #summary of model output
## Family: nbinom2 ( log )
## Formula: dev.seed ~ treat/flowers + (1 | site/plot.treat)
## Zero inflation: ~1
## Data: mert
##
## AIC BIC logLik deviance df.resid
## 783.1 803.7 -383.6 767.1 89
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plot.treat:site (Intercept) 0.02665 0.1633
## site (Intercept) 0.07150 0.2674
## Number of obs: 97, groups: plot.treat:site, 8; site, 4
##
## Overdispersion parameter for nbinom2 family (): 5.6
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.376715 0.242589 9.797 < 2e-16 ***
## treatopen-hand -0.097678 0.255614 -0.382 0.702
## treatopen:flowers 0.024074 0.004401 5.470 4.51e-08 ***
## treatopen-hand:flowers 0.026340 0.004167 6.320 2.61e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.9293 0.4668 -6.275 3.49e-10 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
mert.counts.treat.simulation<-simulateResiduals(fittedModel = mert.treat.glmmtmb) #creating simulation for Mertensia mixed effects model with DHARMa package
plot(mert.counts.treat.simulation) #plotting simulation
testDispersion(mert.counts.treat.simulation) #checking for overdispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.85214, p-value = 0.592
## alternative hypothesis: two.sided
testZeroInflation(mert.counts.treat.simulation) #checking for zero inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0121, p-value = 1
## alternative hypothesis: two.sided
The data are no longer zero-inflated, and we still detected an effect of pollen limitation on the Mertensia seed counts. Mertensia individuals are pollen limited.
Figures
Below is the plot of seed set vs. relative position of blooming in the population for Mertensia individuals. Plot treatment is indicated by individual color.
ggplot(mert.open, aes(x=mert.open$relative.position, y=mert.open$dev.seed, group=mert.open$plot.treat)) +
labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Mertensia Blooming Time on Total Developed Seeds") +
geom_point(aes(color=plot.treat),size=3, alpha = 0.8) +
scale_color_brewer(palette="Dark2") + labs(color="Plot Treatment") +
scale_color_manual(labels= c("control","accelerated"), values = c("control"="darkred","manip"="darkorange")) +
geom_smooth()
ggplot(mert.open, aes(x=mert.open$relative.position, y=mert.open$dev.seed, group=mert.open$early.late)) +
labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Mertensia Blooming Time on Total Developed Seeds") +
theme_classic() +
geom_point(aes(color=early.late),size=3, alpha = 0.8) +
scale_color_brewer(palette="Dark2") + labs(color="Before or After Peak") +
scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange")) +
geom_smooth(method = "glm", aes(color=early.late))
We then plotted the interaction between deviation and early/late instead of relative position in order to fit a linear model.
ggplot(mert.open, aes(x=mert.open$deviation, y=mert.open$dev.seed/mert.open$flowers, group= mert.open$early.late)) +
labs(y="Total Developed Seeds per Flower",x="Days Relative to Population Peak Bloom", title="Effect of Mertensia Deviation from Peak on Total Developed Seeds") + theme_classic() +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=early.late)) +
geom_smooth(method="glm",aes(color=early.late)) +
scale_color_brewer(palette="Dark2") +
labs(color="Before or After Peak") +
scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange"))
Then we plotted the effect of conspecific density and plot treatment on total developed seeds.
ggplot(mert.open, aes(x=mert.open$number.conspecifics, y=mert.open$dev.seed, group= mert.open$plot.treat)) +
labs(y="Total Developed Seeds",x="Conspecific Density", title="Effect of Mertensia Conspecific Density on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=plot.treat)) +
geom_smooth(method="glm",aes(color=plot.treat)) +
scale_color_brewer(palette="Dark2") +
labs(color="Plot Treatment") +
scale_color_manual(labels= c("control","accelerated snowmelt"), values = c("control"="darkgreen","manip"="blue"))
Finally, we plotted all of the data with the pollinator active period. The pollinator active period was the time during which we saw interactions between pollinators and our species. The active period is shaded.
mert$midpoint2<-format(as.Date(mert$midpoint), format="%m-%d")
ggplot(mert, aes(x=mert$midpoint2, y=mert$dev.seed, group= mert$treat)) +
labs(y="Total Developed Seeds",x="Date", title="Effect of Mertensia Blooming Time and Pollinator Activity on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=treat)) +
geom_smooth(aes(color=treat)) +
scale_color_brewer(palette="Dark2") +
labs(color="Open or Hand-pollinated") +
scale_color_manual(labels= c("open","open-hand"), values = c("open"="darkred","open-hand"="darkorange")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
annotate("rect",xmin="06-19", xmax="07-03",ymin=0,ymax=Inf,alpha=0.1)
Delphinium
For Delphinium, we again limited the individuals to open treatment individuals.
delph.open<-delph[(delph$treat=="open"),] #removed hand pollination treatment
delph.nbinom<-fitdist(delph.open$total.seeds, "nbinom") #fitting negative binomial distribution using fitdistplus package
plot(delph.nbinom) #plotting distribution fit
The Delphinium data fits a negative binomial distribution. We used the glmmTMB package again.
Main Model
delph.glmmtmb<-glmmTMB(total.seeds ~ X.flowers:(deviation*early.late) + X.flowers:number.conspecifics/plot.treat + (1|site), family = "nbinom2", data = delph.open) #Delphinium model for negative binomial data with glmmTMB package
#removed treat fixed effect (only open)
summary(delph.glmmtmb) #summary of model output
## Family: nbinom2 ( log )
## Formula:
## total.seeds ~ X.flowers:(deviation * early.late) + X.flowers:number.conspecifics/plot.treat +
## (1 | site)
## Data: delph.open
##
## AIC BIC logLik deviance df.resid
## 403.9 417.6 -193.9 387.9 33
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.02556 0.1599
## Number of obs: 41, groups: site, 4
##
## Overdispersion parameter for nbinom2 family (): 1.37
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) 2.990040 0.369707 8.088
## X.flowers:deviation 0.026479 0.011771 2.249
## X.flowers:early.latelate 0.137212 0.168217 0.816
## X.flowers:number.conspecifics 0.001297 0.001946 0.666
## X.flowers:deviation:early.latelate -0.010662 0.036032 -0.296
## X.flowers:number.conspecifics:plot.treatmanip -0.002071 0.001801 -1.150
## Pr(>|z|)
## (Intercept) 6.09e-16 ***
## X.flowers:deviation 0.0245 *
## X.flowers:early.latelate 0.4147
## X.flowers:number.conspecifics 0.5051
## X.flowers:deviation:early.latelate 0.7673
## X.flowers:number.conspecifics:plot.treatmanip 0.2501
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
delph.counts.simulation<-simulateResiduals(fittedModel = delph.glmmtmb) #creating simulation for Delphinium mixed effects model with DHARMa package
plot(delph.counts.simulation) #plotting simulation
testDispersion(delph.counts.simulation) #checking for overdispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.81708, p-value = 0.496
## alternative hypothesis: two.sided
testZeroInflation(delph.counts.simulation) #checking for zero inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 4.5872, p-value = 0.152
## alternative hypothesis: two.sided
When flower number is taken into account, deviation from the population peak has a significant effect on Delphinium developed seed count. We tested the model assumptions as we did with Mertensia seed count data. The model for the Delphinium count data is not overdispersed or zero-inflated.
To determine the direction of deviation effect, we looked at the plot. When flower number is taken into account, seed set of Delphinium individuals increase with deviation from the population peak.
Pollen Limitation Model
Again, we tested the pollen limitation in a separate mixed effects model.
#delph.treat.glmmtmb<-glmmTMB(total.seeds ~ treat/X.flowers + (1|site/plot.treat), family = "nbinom2", data = delph) #Mertensia model for negative binomial data with glmmTMB package
#summary(delph.treat.glmmtmb) #summary of model output
#delph.counts.treat.simulation<-simulateResiduals(fittedModel = delph.treat.glmmtmb) #creating simulation for Delphinium mixed effects model with DHARMa package
#plot(delph.counts.treat.simulation) #plotting simulation
#testDispersion(delph.counts.treat.simulation) #checking for overdispersion
#testZeroInflation(delph.counts.treat.simulation) #checking for zero inflation
We detected an effect of pollen limitation on the total developed seeds of Delphinium individuals. We again checked for overdispersion and zero-inflation in the pollen limitation data. The Delphinum pollen limitation model is overdispersed and zero-inflated. We adjusted our model to include the zero-inflation formula, and re-tested for overdispersion and zero-inflation.
delph.treat.glmmtmb<-glmmTMB(total.seeds ~ treat/X.flowers + (1|site/plot.treat), ziformula=~1, family = "nbinom2", data = delph) #Mertensia model for negative binomial data with glmmTMB package
summary(delph.treat.glmmtmb) #summary of model output
## Family: nbinom2 ( log )
## Formula: total.seeds ~ treat/X.flowers + (1 | site/plot.treat)
## Zero inflation: ~1
## Data: delph
##
## AIC BIC logLik deviance df.resid
## 937.0 957.5 -460.5 921.0 88
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plot.treat:site (Intercept) 0.003478 0.05897
## site (Intercept) 0.071347 0.26711
## Number of obs: 96, groups: plot.treat:site, 7; site, 4
##
## Overdispersion parameter for nbinom2 family (): 1.97
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 2.54661 0.36803 6.920 4.53e-12 ***
## treatopen-hand 0.55500 0.47352 1.172 0.24117
## treatopen:X.flowers 0.29208 0.06842 4.269 1.97e-05 ***
## treatopen-hand:X.flowers 0.18467 0.06693 2.759 0.00579 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.2789 0.5892 -5.565 2.63e-08 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
delph.counts.treat.simulation<-simulateResiduals(fittedModel = delph.treat.glmmtmb) #creating simulation for Delphinium mixed effects model with DHARMa package
plot(delph.counts.treat.simulation) #plotting simulation
testDispersion(delph.counts.treat.simulation) #checking for overdispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.73169, p-value = 0.096
## alternative hypothesis: two.sided
testZeroInflation(delph.counts.treat.simulation) #checking for zero inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0604, p-value = 1
## alternative hypothesis: two.sided
The new model is not overdispersed or zero-inflated. We still see an effect of pollen limitation on the Delphinium count data. Delphinium individuals are pollen limited.
Figures
Below are the plots for the Delphinium individuals.
ggplot(delph.open, aes(x=delph.open$relative.position, y=delph.open$total.seeds, group=plot.treat)) +
labs(title="Effect of Delphinium Blooming Time on Total Developed Seeds", y="Total Developed Seeds",x="Days Relative to Population Peak Bloom") +
geom_point(aes(color=plot.treat),size=3, alpha = 0.8) +
scale_color_brewer(palette="Dark2") + labs(color="Plot Treatment") +
scale_color_manual(labels= c("control","accelerated"), values = c("control"="darkred","manip"="darkorange")) +
geom_smooth()
ggplot(delph.open, aes(x=delph.open$relative.position, y=delph.open$total.seeds, group=delph.open$early.late)) +
labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Delphinium Blooming Time on Total Developed Seeds") +
theme_classic() +
geom_point(aes(color=early.late),size=3, alpha = 0.8) +
scale_color_brewer(palette="Dark2") + labs(color="Before or After Peak") +
scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange")) +
geom_smooth(method = "glm", aes(color=early.late))
ggplot(delph.open, aes(x=delph.open$deviation, y=(delph.open$total.seeds/delph.open$X.flowers), group= delph.open$early.late)) +
labs(y="Total Developed Seeds per Flower",x="Days Relative to Population Peak Bloom", title="Effect of Delphinium Deviation from Peak on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=early.late)) +
geom_smooth(method="glm",aes(color=early.late)) +
scale_color_brewer(palette="Dark2") +
labs(color="Before or After Peak") +
scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange"))
ggplot(delph.open, aes(x=delph.open$number.conspecifics, y=delph.open$total.seeds, group= delph.open$plot.treat)) +
labs(y="Total Developed Seeds",x="Conspecific Density", title="Effect of Delphinium Conspecific Density on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=plot.treat)) +
geom_smooth(method="glm",aes(color=plot.treat)) +
scale_color_brewer(palette="Dark2") +
labs(color="Plot Treatment") +
scale_color_manual(labels= c("control","accelerated snowmelt"), values = c("control"="darkgreen","manip"="blue"))
Finally, we plotted the pollinator active period with the total number of developed seeds and midpoint of individuals.
delph$midpoint2<-format(as.Date(delph$midpoint), format="%m-%d")
ggplot(delph, aes(x=delph$midpoint2, y=delph$total.seeds, group= delph$treat)) +
labs(y="Total Developed Seeds",x="Date", title="Effect of Delphinium Blooming Time and Pollinator Activity on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=treat)) +
geom_smooth(aes(color=treat)) +
scale_color_brewer(palette="Dark2") +
labs(color="Open or Hand-pollinated") +
scale_color_manual(labels= c("open","open-hand"), values = c("open"="darkred","open-hand"="darkorange")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
annotate("rect",xmin="06-27", xmax="07-18",ymin=0,ymax=Inf,alpha=0.1)
Potentilla
We began by limiting the data to the open treatment Potentilla individuals.
pot.open<-pot[(pot$treat=="open"),] #removed hand pollination treatment
pot.nbinom<-fitdist(pot.open$total.seeds, "nbinom") #fitting negative binomial distribution using fitdistplus package
plot(pot.nbinom) #plotting distribution fit
The Potentilla seed data fits a negative binomial distribution. We used the glmmTMB package again.
Main Model
#pot.glmmtmb<-glmmTMB(total.seeds ~ X.flowers:(deviation * early.late) + X.flowers:number.conspecifics/plot.treat + (1|site), family = "nbinom2", data = pot.open) #Potentilla model for negative binomial data with glmmTMB package
#summary(pot.glmmtmb) #summary of model output
#pot.counts.simulation<-simulateResiduals(fittedModel = pot.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
#plot(pot.counts.simulation) #plotting simulation
#testDispersion(pot.counts.simulation) #checking for overdispersion
#testZeroInflation(pot.counts.simulation) #checking for zero inflation
The Potentilla data are not overdispersed, but they are zero-inflated. To fix this, we had to include zero inflation in our model and again check for zero-inflation.
pot.glmmtmb<-glmmTMB(total.seeds ~ X.flowers:(deviation * early.late) + X.flowers:number.conspecifics/plot.treat + (1|site), ziformula=~1, family = "nbinom2", data = pot.open) #new Potentilla model for negative binomial data with glmmTMB package
summary(pot.glmmtmb) #summary of model output
## Family: nbinom2 ( log )
## Formula:
## total.seeds ~ X.flowers:(deviation * early.late) + X.flowers:number.conspecifics/plot.treat +
## (1 | site)
## Zero inflation: ~1
## Data: pot.open
##
## AIC BIC logLik deviance df.resid
## 280.7 291.4 -131.4 262.7 15
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## site (Intercept) 3.537e-10 1.881e-05
## Number of obs: 24, groups: site, 4
##
## Overdispersion parameter for nbinom2 family (): 2.95
##
## Conditional model:
## Estimate Std. Error z value
## (Intercept) 3.657e+00 4.095e-01 8.930
## X.flowers:deviation 1.220e-02 4.421e-03 2.760
## X.flowers:early.latelate 7.605e-02 2.733e-02 2.783
## X.flowers:number.conspecifics 1.631e-05 3.236e-04 0.050
## X.flowers:deviation:early.latelate -1.423e-02 4.629e-03 -3.074
## X.flowers:number.conspecifics:plot.treatmanip 4.595e-04 3.227e-04 1.424
## Pr(>|z|)
## (Intercept) < 2e-16 ***
## X.flowers:deviation 0.00578 **
## X.flowers:early.latelate 0.00539 **
## X.flowers:number.conspecifics 0.95981
## X.flowers:deviation:early.latelate 0.00211 **
## X.flowers:number.conspecifics:plot.treatmanip 0.15444
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -3.138 1.024 -3.065 0.00218 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pot.counts.simulation<-simulateResiduals(fittedModel = pot.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
plot(pot.counts.simulation) #plotting simulation
testDispersion(pot.counts.simulation) #checking for overdispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.8079, p-value = 0.52
## alternative hypothesis: two.sided
testZeroInflation(pot.counts.simulation) #checking for zero inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0638, p-value = 1
## alternative hypothesis: two.sided
When we incorporated zero-inflation into the model, we found a significant effect of deviation from the population peak, blooming before vs after the peak, and the interaction between the two. The direction of the effect, however, must be determined from the plots (below).
After looking at the plots, we see that when we account for the number of flowers, deviation has a greater negative effect on the late blooming Potentilla individuals.
Pollen Limitation Model
Below is the mixed effects model for the pollen limitation treatments.
#pot.treat.glmmtmb<-glmmTMB(total.seeds ~ treat/X.flowers + (1|site/plot.treat), family = "nbinom2", data = pot) #Potentilla model for negative binomial data with glmmadmb package
#summary(pot.treat.glmmtmb) #summary of model output
#pot.counts.treat.simulation<-simulateResiduals(fittedModel = pot.treat.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
#plot(pot.counts.treat.simulation) #plotting simulation
#testDispersion(pot.counts.treat.simulation) #checking for overdispersion
#testZeroInflation(pot.counts.treat.simulation) #checking for zero inflation
We detected a significant effect of the pollen limitation treatment on the total number of developed seeds for Potentilla individuals, but the Potentilla pollen limitation data were overdispersed and zero-inflated. We addressed the zero-inflation issue by again including the zero-inflation formula into our models.
pot.treat.glmmtmb<-glmmTMB(total.seeds ~ treat/X.flowers + (1|site/plot.treat), ziformula=~1, family="nbinom2",data = pot) #Potentilla model for negative binomial data with glmmadmb package
summary(pot.treat.glmmtmb) #summary of model output
## Family: nbinom2 ( log )
## Formula: total.seeds ~ treat/X.flowers + (1 | site/plot.treat)
## Zero inflation: ~1
## Data: pot
##
## AIC BIC logLik deviance df.resid
## 502.1 516.4 -243.0 486.1 36
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev.
## plot.treat:site (Intercept) 6.676e-09 8.171e-05
## site (Intercept) 1.193e-02 1.092e-01
## Number of obs: 44, groups: plot.treat:site, 6; site, 4
##
## Overdispersion parameter for nbinom2 family (): 1.86
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 3.87623 0.56142 6.904 5.04e-12 ***
## treatopen-hand 0.61864 0.68954 0.897 0.3696
## treatopen:X.flowers 0.05534 0.03138 1.764 0.0778 .
## treatopen-hand:X.flowers 0.01014 0.02816 0.360 0.7189
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -2.6271 0.6051 -4.342 1.41e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
pot.counts.treat.simulation<-simulateResiduals(fittedModel = pot.treat.glmmtmb) #creating simulation for Potentilla mixed effects model with DHARMa package
plot(pot.counts.treat.simulation) #plotting simulation
testDispersion(pot.counts.treat.simulation) #checking for overdispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.80175, p-value = 0.336
## alternative hypothesis: two.sided
testZeroInflation(pot.counts.treat.simulation) #checking for zero inflation
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0204, p-value = 1
## alternative hypothesis: two.sided
The data are no longer zero-inflated or overdispersed, and we still detected a marginally significant effect of pollen limitation in the Potentilla individuals. Potentilla individuals are slightly pollen limited.
Figures
Below are the plots for Potentilla individuals.
ggplot(pot.open, aes(x=pot.open$relative.position, y=pot.open$total.seeds, group=plot.treat)) +
labs(title="Effect of Potentilla Blooming Time on Total Developed Seeds", y="Total Developed Seeds",x="Days Relative to Population Peak Bloom") +
geom_point(aes(color=plot.treat),size=3, alpha = 0.8) +
scale_color_brewer(palette="Dark2") + labs(color="Plot Treatment") +
scale_color_manual(labels= c("control","accelerated"), values = c("control"="darkred","manip"="darkorange")) +
geom_smooth()
ggplot(pot.open, aes(x=pot.open$relative.position, y=pot.open$total.seeds, group=pot.open$early.late)) +
labs(y="Total Developed Seeds",x="Days Relative to Population Peak Bloom", title="Effect of Potentilla Blooming Time on Total Developed Seeds") +
theme_classic() +
geom_point(aes(color=early.late),size=3, alpha = 0.8) +
scale_color_brewer(palette="Dark2") + labs(color="Before or After Peak") +
scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange")) +
geom_smooth(method = "glm", aes(color=early.late))
ggplot(pot.open, aes(x=pot.open$deviation, y=pot.open$total.seeds/pot.open$X.flowers, group= pot.open$early.late)) +
labs(y="Total Developed Seeds per Flower",x="Days Relative to Population Peak Bloom",title="Effect of Potentilla Deviation from Peak on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=early.late)) +
geom_smooth(method="glm",aes(color=early.late)) +
scale_color_brewer(palette="Dark2") +
labs(color="Before or After Peak") +
scale_color_manual(labels= c("before peak","after peak"), values = c("early"="darkred","late"="darkorange"))
ggplot(pot.open, aes(x=pot.open$number.conspecifics, y=pot.open$total.seeds, group= pot.open$plot.treat)) +
labs(y="Total Developed Seeds",x="Conspecific Density", title="Effect of Potentilla Conspecific Density on Total Developed Seeds") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=plot.treat)) +
geom_smooth(method="glm",aes(color=plot.treat)) +
scale_color_brewer(palette="Dark2") +
labs(color="Plot Treatment") +
scale_color_manual(labels= c("control","accelerated snowmelt"), values = c("control"="darkgreen","manip"="blue"))
Again we plotted the pollinator activity period.
ggplot(pot,aes(x=as.Date.factor(pot$midpoint), y=pot$total.seeds, group= pot$treat)) +
labs(y="Total Developed Seeds",x="Date", title="Pollinator Activity Period and Total Developed Seeds in Potentilla") +
theme_classic() +
geom_point(size=3, alpha = 0.8, aes(color=treat)) +
geom_smooth(aes(color=treat)) +
scale_color_brewer(palette="Dark2") +
labs(color="Open or Hand-pollinated") +
scale_color_manual(labels= c("open","open-hand"), values = c("open"="darkred","open-hand"="darkorange")) +
theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
annotate("rect",xmin=as.Date.character("2019-07-05"), xmax=as.Date.character("2019-08-24"),ymin=0,ymax=Inf,alpha=0.1) +
scale_x_date(date_breaks = "3 days", date_labels = "%m/%d")
Mixed Effects Models with Seeds per Flower
Mertensia
mert.open$non_zero <- ifelse(mert.open$seeds.per.flower > 0, 1, 0) # creating zero/non-zero column
# also running into errors when evaluating this first line. Error: "Response is constant." Could be from not enough 0 entries?
mert.hurdle<-glmer(non_zero ~ early.late*deviation + number.conspecifics/plot.treat + rate + (1|site), data=mert.open,family = binomial) # model with new column
mert.open.nonzero = subset(mert.open,non_zero==1)
mert.gamma<-glmer(seeds.per.flower ~ early.late*deviation + number.conspecifics/plot.treat + rate + (1|site), family=Gamma, data = mert.open.nonzero, na.action = "na.fail") # model without zeros
summary(mert.hurdle) # viewing model results
summary(mert.gamma)
None of the variables affected the presence or absence of seeds for Mertensia individuals. When evaluating the individuals that produced seeds, we detect an effect of conspecific density and soil moisture decline on the seeds per flower.
Delphinium
delph.open$non_zero <- ifelse(delph.open$seeds.per.flower > 0, 1, 0) # creating zero/non-zero column
delph.hurdle<-glmer(non_zero ~ early.late*deviation + number.conspecifics/plot.treat + rate + (1|site), data=delph.open,family = binomial) # model with new column
delph.open.nonzero = subset(delph.open,non_zero==1)
delph.gamma<-glmer(seeds.per.flower ~ early.late*deviation + number.conspecifics/plot.treat + rate + (1|site), family=Gamma, data = delph.open.nonzero) # model without zeros
summary(delph.hurdle) # viewing model results
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: binomial ( logit )
## Formula: non_zero ~ early.late * deviation + number.conspecifics/plot.treat +
## rate + (1 | site)
## Data: delph.open
##
## AIC BIC logLik deviance df.resid
## 20.2 28.9 -2.1 4.2 14
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -0.001014 0.000000 0.000042 0.000205 0.041321
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 31179 176.6
## Number of obs: 22, groups: site, 4
##
## Fixed effects:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 4.885e+01 4.314e+03 0.011 0.991
## early.latelate -7.407e+04 7.925e+07 -0.001 0.999
## deviation -4.492e+00 7.054e+02 -0.006 0.995
## number.conspecifics 7.282e-01 2.736e+02 0.003 0.998
## rate 8.626e+00 2.694e+02 0.032 0.974
## early.latelate:deviation 7.521e+04 1.886e+07 0.004 0.997
## number.conspecifics:plot.treatmanip -8.663e-01 2.683e+02 -0.003 0.997
##
## Correlation of Fixed Effects:
## (Intr) erly.l devitn nmbr.c rate erly.:
## early.latlt 0.000
## deviation -0.766 0.000
## nmbr.cnspcf -0.871 0.000 0.523
## rate 0.349 0.000 0.304 -0.458
## erly.ltlt:d 0.000 -0.872 0.000 0.000 0.000
## nmbr.cnsp:. 0.898 0.000 -0.628 -0.990 0.361 0.000
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## unable to evaluate scaled gradient
## Model failed to converge: degenerate Hessian with 2 negative eigenvalues
summary(delph.gamma)
## Generalized linear mixed model fit by maximum likelihood (Laplace
## Approximation) [glmerMod]
## Family: Gamma ( inverse )
## Formula:
## seeds.per.flower ~ early.late * deviation + number.conspecifics/plot.treat +
## rate + (1 | site)
## Data: delph.open.nonzero
##
## AIC BIC logLik deviance df.resid
## 151.0 160.0 -66.5 133.0 11
##
## Scaled residuals:
## Min 1Q Median 3Q Max
## -1.96294 -0.26010 0.01632 0.42403 2.39969
##
## Random effects:
## Groups Name Variance Std.Dev.
## site (Intercept) 0.000 0.000
## Residual 0.243 0.493
## Number of obs: 20, groups: site, 4
##
## Fixed effects:
## Estimate Std. Error t value Pr(>|z|)
## (Intercept) 0.1205735 0.0668481 1.804 0.07128 .
## early.latelate -0.0046279 0.0911840 -0.051 0.95952
## deviation 0.0021853 0.0111990 0.195 0.84529
## number.conspecifics 0.0002483 0.0006457 0.385 0.70061
## rate 0.0139848 0.0049891 2.803 0.00506 **
## early.latelate:deviation -0.0097170 0.0176989 -0.549 0.58299
## number.conspecifics:plot.treatmanip -0.0005570 0.0006800 -0.819 0.41272
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Correlation of Fixed Effects:
## (Intr) erly.l devitn nmbr.c rate erly.:
## early.latlt -0.719
## deviation -0.833 0.674
## nmbr.cnspcf -0.340 0.187 0.134
## rate -0.086 0.209 0.470 0.117
## erly.ltlt:d 0.513 -0.904 -0.650 -0.034 -0.343
## nmbr.cnsp:. 0.116 -0.102 -0.373 -0.525 -0.383 0.226
## optimizer (Nelder_Mead) convergence code: 0 (OK)
## boundary (singular) fit: see ?isSingular
As with Mertensia, none of the variables affected the ability for individuals to produce seeds. In terms of the Delphinium individuals that did produce seeds, the soil moisture decline has a significant effect on the seeds per flower in Delphinium.
Potentilla
pot.open$non_zero <- ifelse(pot.open$seeds.per.flower > 0, 1, 0) # creating zero/non-zero column
# had issues running these ones because so few open treatment individuals were 0 in the zero/non-zero breakdown, got error message and I'm still working on it
pot.hurdle<-glmer(non_zero ~ early.late*deviation + number.conspecifics/plot.treat + rate + (1|site), data=pot.open,family = binomial) # model with new column
pot.open.nonzero = subset(pot.open,non_zero==1)
pot.gamma<-glmer(seeds.per.flower ~ early.late*deviation + number.conspecifics/plot.treat + rate + (1|site), family=Gamma, data = pot.open.nonzero) # model without zeros
summary(pot.hurdle) # viewing model results
summary(pot.gamma)
